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ABSTRACT 

We present measurements of the clustering of hot and cold patches in the mi- 
crowave background sky as measured from the Wilkinson Microwave ^nisotropy 
Probe (WMAP) five-year data. These measurements are compared with theoretical 
predictions which assume that the cosmological signal obeys Gaussian statistics. We 
find significant differences from the simplest Gaussian-based prediction. However, the 
measurements are sensitive to the fact that the noise is spatially inhomogeneous (e.g., 
because different parts of the sky were observed for different lengths of time) . We show 
how to account for this spatial inhomogeneity when making predictions. Differences 
from the Gaussian-based expectation remain even after this more careful accounting 
of the noise. In particular, we note that hot and cold pixels cluster differently within 
the same temperature thresholds at few-degree scales. While these findings may in- 
dicate primordial non-Gaussianity, we discuss other plausible explanations for these 
discrepancies. In addition, we find some deviations from Gaussianity at sub-degree 
scales, especially in the W band, whose origin may be associated with extragalactic 
dust emission. 

Key words: methods: statistical - methods: data analysis - cosmology: cosmic mi- 
crowave background, correlations — cosmology: observations. 



1 INTRODUCTION 

Non-Gaussianity of the Cosmic Microwave Background 
(CMB) is progressively becoming a crucial probe of the in- 
flationary dynamics, since models beyond single field slow- 
roll inflation predict large distinct non-Gaussianities above 
the current limit of measurement (Lyth et al. 2003, 2006; 
Bartolo et al. 2004; Dvafi et al. 2004; Chen 2005; Seery & 
Lidsey 2005; Chen & Szapudi 2006; Seery & Hidalgo 2006; 
Ling & Wu 2008; Senatore et al. 2008, 2009; Chen et al. 
2009). Moreover, up-to-date a number of plausible alterna- 
tive early universe models also predict skewed primordial 
fluctuations, and are in principle distinguishable from other 
scenarios through the shape dependence of high-order cor- 
relation functions (Alishahiha et al. 2004; Arkani-Hamed et 
al. 2004; CremineUi & Senatore 2007; Brandenberger 2008; 
Buchbinder, Khoury & Ovrut 2008; Lehners & Steinhardt 
2008; McAllister & Silverstein 2008). 

An even small degree of primordial non-Gaussianity 
would indicate a quite different structure formation scenario 
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from the concordance cosmological model in which density 
perturbations are assumed to be a Gaussian random fleld, 
and alter significantly the statistics of voids (Kamionkowski 
et al. 2008; Song & Lee 2008, for example). Imprints of 
primordial non-Gaussianity can manifest in the statistical 
properties of the Lyman-a forest QSO spectra at intermedi- 
ate redshift (Viel et al. 2008), in the large-scale distribution 
of neutral hydrogen (Pillepich et al. 2007), in the reioniza- 
tion history of the universe (Crociani et al. 2008), and in 
the abundance, clustering and biasing of dark matter halos 
(Kang et al. 2007; Carbone et al. 2008; Dalai et al. 2008; 
Desjacques et al. 2008; Grossi et al. 2008; LoVerde et al. 
2008; Matarrese & Verde 2008; McDonald & Seljak 2008; 
Pillepich et al. 2008; Seljak 2008; Slosar et al. 2008; Jeong & 
Komatsu 2009). Topology of the large scale structure is also 
very sensitive to the initial conditions of the matter density 
field, and potentially represents a strong and independent 
test of deviations from the Gaussian hypothesis (Park et al. 
1998, 2001, 2005; Gott et al. 2008a,b; Hikage et al. 2008b). 

Hence, there is significant interest in quantifying the 
Gaussianity of the cosmic background radiation, and a 
multitude of non-Gaussian estimators have been applied 
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so far. In particular, data from the lyilkinson Microwave 
Anisotropy Probe (WMAP) mission has been used to con- 
strain primordial non-Gaussianity. Recent measurements of 
the one- and three-point probabihty functions, the bispec- 
trum and trispectrum, the genus statistic and the other 
Minkowski functionals (Cabella et al. 2005; Hinshaw et al. 
2007, 2008; CremineUi et al. 2007; Gott et al. 2007; Spergel 
et al. 2007; Curto et al. 2008; Hikage et al. 2008a; Komatsu 
et al. 2008; Vielva & Sanz 2008; Pietrobon et al. 2009a) from 
the first three or five years of WMAP data are all consis- 
tent with Gaussianity. However, a considerable number of 
anomalies like asymmetries, intrinsic alignments in the data 
or the presence of unusual cold spots have also been reported 
(Cruz et al. 2006; Tojeiro et al. 2006; Chiang et al. 2007; 
Naselsky et al. 2007; Rath et al. 2007, 2009; Vielva et al. 
2007; Hansen et al. 2008; Dickinson et al. 2009; Diego et al. 
2009; Li et al. 2009; Pietrobon et al. 2009b; Rossmanith et 
al. 2009), as well as some claims of non-Gaussianity (Jeong 
& Smoot 2007; Yadav & Wandelt 2008), which, if confirmed, 
would rule out a large class of inflationary models. 

In this paper we investigate the Gaussian hypothesis us- 
ing the clustering statistics of pixels that lie above or below a 
threshold. For a Gaussian field, the scale dependence of this 
statistic, and its dependence on threshold, have been pre- 
dicted (Jensen & Szalay 1986). It turns out that a careful ac- 
counting for the effects of noise is crucial to performing this 
test. Kashlinsky et al. (2001) and Hernandez-Monteagudo 
et al. (2004) presented an analysis of the problem when the 
noise is spatially homogeneous. In the case of WMAP, the 
noise is spatially inhomogeneous; the present work shows 
how to incorporate this complexity into the analysis. 

The layout of the paper is as follows. In Section [2] we 
outline the theoretical framework for the two-point statistics 
above a temperature threshold, and show how to incorpo- 
rate the complexity of inhomogeneous noise into the two- 
point formalism. Some details on errorbar calculations are 
left in Appendix A. In Section Owe test our theory against 
WMAP five-year data. We study and quantify the inhomo- 
geneous noise properties and show that the abundance and 
clustering of pixels (both scale and dependence on thresh- 
old temperature) appears to be inconsistent with the Gaus- 
sian prediction if one ignores the fact that the noise is spa- 
tially inhomogeneous. We compare the number density of 
pixels above threshold and the one- and two-point statistics 
measurements with the theory, allowing for inhomogeneous 
noise. While this improves the agreement between measure- 
ments and predictions, discrepancies remain. Although this 
may be an indication of primordial non-Gaussianity, we dis- 
cuss a variety of other possibilities, among which the effect 
of smoothing the map. We inspect all the WMAP five-year 
channels, but we present in the main text results only for 
the Wl differencing assembly (DA), and leave in Appendix 
B results for other DA's. A final section summarizes our 
findings, and highlights ongoing and future studies. 



2 THEORETICAL MODEL 
2.1 Basic Notation 

Denote the observed value in a pixel by _D = T — (T) = 5T, 
which is the sum of the true signal s plus noise n, both 



of which have mean zero. We consider a model in which 
the signal is homogeneous and may have spatial correlations 
whereas the noise, which is independent of the signal, may 
be inhomogeneous and have spatial correlations. By this we 
mean that the rms value of the noise (Jn may fiuctuate from 
pixel to pixel and these fiuctuations may be correlated, and 
that the actual value of the noise in one pixel could depend 
of that in another. Let p{D) denote the observed one-point 
distribution of D, G(s) the distribution of s, p(o"n) the dis- 
tribution of the rms value of the noise in a pixel, and (;(n|(Tn) 
the distribution of the noise when the rms value of the noise 
is (Tn. Note that g{n) — J g{n\an)p{o'n)dan- Our convention 
is to use capital letters for average quantities, and lower case 
letters for actual (varying) quantities. Later on in the paper 
we will also drop the (understood) subscript n of the rms of 
the noise, for clarity of notation. 



2.2 One-Point Observed Distribution 

Since the distribution of the noise g{n) is independent of the 
signal s, 



p{D) 



ds G(s) / An g{n)5ji{s + n = D) 



ds G{s) J '^^ J 'i'^n 5'("-kn) P(o-n) 

X 5r)(s + n = D) 

dan p{o-n) / ds G(s) g{D - s\an) 



d(7n P(crn) p{D\an) 

where So is the Dirac delta and 
p{D\an) = y ds G(s) g(D - s\an). 
The variance of D is; 
al = (D^) = / dDp{D)D^ 



(1) 



(2) 



danp{an) / dsG(s) / dn g{n\an) {s + n) 



al+ danp{an) {n^\an} 



2 , 2 



(3) 



where we have used as to denote the rms value of the signal, 
and a^ to denote the variance of the noise upon averaging 
over all pixels. It is also straightforward to show that (D) — 
0, as expected. 

The fraction of pixels above some threshold Dt is 

/•oo p poo 

/(A) = / dDp{D) = / danp{an) / dDp{D\an) 

Jot J JDt 



danp{an)f{Dt\an). (4) 
If {D\Dt} denotes the mean value of D in such pixels, then 

/• oo 

f{D,){D\D,) = / dDp{D)D 

JDt 
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/poci 
dfTnp(o-„) j dDp{D\ar,)D 

= j danPian)f{Dt\an){D\Dt,an}. (5) 



2.3 Two-Point Observed Distribution 

Two point statistics may be computed similarly. For two 
pixels separated by the angular distance 6, 



p{Di,D2\e) = J dsi J ds2 G(si,S2\e) J dm J dn2 
X g(ni,n2) (5d(si + ni = Di) (5d(,S2 + n2 = D2) 
= j dsi j ds2 G{si,S2\d) j dni j dn2 j dai 

X j diT2 5(ni,n2|fTi,a-2) p(fTi,cT2|^) i5d(si + ni = Di) 
X (5d(s2 + n2 = D2) 

= j dai j da2p{(Ti,a2\0) j dsi j ds2G{si,S2\0) 

X - si |cti) g{D2 - S2\(J2) 

= j dai j da2p{cTi,cT2\e)p{Di,D2\ai,a2,e) (6) 

where 

p{Di,D2\au(72,9) = J d.si J ds2G{si,S2\e) 

X g{Di - si\cn) g{D2 - 52^2). (7) 
Notice that if we integrate over all pixels then 



{DiD2\e) = y dDi y dD2DiD2piDi,D2\e) 
= J dsi J ds2G{si,S2\0) S1S2 
+ / dni / dn2 g(ni,n2\9) nin2 



= {siS2\e} + {nin2\e) 



(8) 



where 

{siS2\e)= j dsi j ds2G{si,S2\9)siS2=Cs{0) (9) 

and 

(nin2\9) = J dni J dn2 g{ni,n2\9) nin2 

dni / dn2 / dai / da2 p{ai,a2\9) 

X g{ni,n2\ai,a2,9) nin2 

dcTi / dcr2 p(cri, o-2|6') / dni I dn2 
X g{ni,n2\ai,a2,9) nin2 

oo poo 



Jo Jo 



da2 p{ai,a2\9) {nin2\ai,a2,9) 
= CNi9). (10) 



This is because the signal and noise are uncorrclated, and so 
{sn\9} — for all 6 > 0. Clearly, if the noise is not spatially 
correlated, then {nin2\9} = and (DiD2\9} = Gs{9). 

If the rms values of the noise are perfectly correlated, 
meaning tJi = CT2, then it is convenient to write p{ai,a2\9) = 
p{ai) p{a2\ai , 9) and replace p{a2\ai,9) with a delta func- 
tion centered on ai = an- One then averages over the dis- 
tribution of an, namely: 



p{Di,D2\9) = y^dui J da2p{ai) do{ai = a2 = a^ 



X p{Di,D2\ai,a2,9) 

J dan p{an) p{Di,D2\an,an,9). 



(11) 



If, in addition, one replaces p{an) with a delta function cen- 
tered on (Tn, then p{Di,D2\9) = p{Di,D2\aj^,ai<i,9), and 
the homogeneous case is recovered. The noise distributions 
p{an) andp((Ti,(T2|6') are measured directly from the WMAP 
five-year data. 



2.4 Two-Point Function above Threshold and 
Inhomogeneous Noise 

The two-point unweighted correlation function ^{9) esti- 
mates the excess of probability of finding a pair at a dis- 
tance 9, compared to a random catalog (e.g. Jensen & Sza- 
lay, 1986). Similarly, one can define the two-point weighted 
correlation function W{9), which differs from ^(6) only be- 
cause each member of the pair is now weighted by some (av- 
erage) mark. In this study, the mark is given by ST/{5T) = 
D/(D). Two-point correlation functions of regions above 
some threshold are readily obtained by integration above 
the threshold level. 

Within the formalism previously outlined, the two point 
correlation function of pixels above some threshold Dt is 



l + 6(e) 



dDi 



dD: 



p{Di,D2\9) 



and the associated weighted function is 

D1D2 p{Di,D2\9) 



1 + Wt(6l) 



POO poo 
JOt JDt 



dD2 



(12) 



(13) 



where p{Di, D2\9), {D\Dt) and /(-Dt) have been previously 
defined. 

To gain intuition, it is helpful to define 



l + ^t{9\ai,a2) -- 
and 

1 + Wt(9\ai,a2) 



in terms of which 



J Di J L 



p{Di,D2\ai,a2,9) 
ui ' ^ /(AkO/(Ak2) 



/ D D- 

{D\D.,ai){D\D.,a2) 

p{Di,D2\ai,a2,9) 



/(Aki)/(Ak2)' 



(15) 



1 + ^9) 



dai / da2p{ai,a2\9) 



f{Dt\ai) f{D,\a2) 



/(A)= 



X 1-1- 6(^^1,0-2) 



(16) 
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and 

l + Wt{6) 



dcTl j A(T2P{01,02 

f{DtW2){D\Dt,(J2) 
f{D,){D\D,) 



f{Dt\ai){D\Dt,(7i) 



f{D,) {D\D,) 



1 + W,{e\a^,a2] 



(17) 



This shows that and Wt can be thought of as weighted 
averages over the values at fixed cti and 02- This re- 
writing shows clearly that the key quantity of interest are 
p{D\,D2\(Ji,a2,0), which is the convolution of the signal 
and noise distributions, and p(_D|crn). This rewriting also al- 
lows one to study two limiting cases. If the rms values of the 
noise are perfectly correlated, meaning g\ = a2 = crn, then: 



dDi / dL»2 

Dt J Dt 

p{Dl,D2\an,(Jn,0) 



and 

1 + Wti9) = 



/(A)2 

dcrnp(crn) 



(18) 



D1D2 



p{Dl, D2\an, CTn, 9) 



(19) 



If, in addition, one replaces p(crn) with a delta function cen- 
tered on ctn, then 



l + 6(e) = 
and 

i + Wt{e) 



dDi / dD2 



p(Dl,-D2|(7-N,0-N,^) 

7(A? 



(20) 



DID2 



p{Di, D2\a]^,(7N,d) 



(21) 



We are particularly interested in the case where the 
signal G{si,S2\ff) is bivariate Gaussian with {s\) = (sj) = 
a|, and {siS2\d) = Cs{6) and the noise g{n\a^) is Gaussian 
with rms (Tn. Then 



p{Di,D2\(Jl,(J2,i 

1 



1 



(22) 



27ra2 y; 



exp 



Q1Q2 — Wq 

2 I 



C12DI + aiDl - 2w0 D1D2 
20-2 (aia2 - ™|) 



D l."i"2 UJgJ 

2 „ /_2 I _2\ ;_2 



with Qi = ((Tg -I- o-i)/cri3, Q2 = (cTg -I- cr2)/cri3, and 



■We 



(23) 



+ VK.^^^P M^;-°°"> po (cos 9) , (24) 



47r 



Cn(6') = V i?i±il wr""'"" P°(cos6i), 
47r 



(25) 



where C is the covariance matrix of the temperature 
field, Ci' the power spectrum of the noise map, ^wmap 



the WMAP window function and W 



the additional 



smoothing due to finite pixel size, mask influence and an 
optional additional Gaussian beam smoothing. 



If the noise is spatially uncorrelated, then clearly 
Cn(^) ~ and therefore we = Cs(6) / a^. 

In the approximation where cri =02, rms noise varies 
spatially on scales much larger than those of interest, then 
Qi = a2. The "standard" approximation, rms noise inde- 
pendent of position, has ai = 02 = 1 and we recover the 
Hernandez-Monteagudo et al. (2004) formula. 

Uncertainties in the correlation functions above thresh- 
old are estimated from the optimal variance limit, containing 
cosmic variance, instrumental noise (Knox 1995; Hernandez- 
Monteagudo et al. 2004), and finite binsize effects. See Ap- 
pendix A for more details. 



3 ANALYSIS OF WMAP 5- YEAR DATA 
3.1 Dataset and Pipeline 

The data from the five years of the WMAP mission are avail- 
able online at Legacy Archive for Microwave Background 
Dylta (LAMBDA) website (http://lambda.gsfc.riasa.gov/), 
NASA's CMB Thematic Data Center. For the purposes 
of this study, we use the WMAP five-year full resolution 
"foreground-reduced" coadded sky maps (produced by per- 
forming a weighted pixel-by-pixel mean of the five single year 
maps, see Hinshaw et al. 2008) for the W (the 94 GHz chan- 
nel), V (the 61 GHz channel) and Q (the 41 GHz channel) 
differencing assemblies. We show results for the Wl channel 
in the main text, and leave in Appendix B major results for 
all the other individual DA's. The data consist of four fields 
for each pixel: 

(i) the thermodynamic temperature in mK, 

(ii) the Q Polarization temperature in mK, 

(iii) the U Polarization temperature in mK, and 

(iv) the effective number of observations, A'obs. 

The maps are provided in the HEALPix scheme (Gorski 
et al. 1999) at a resolution of A^side = 512, giving a to- 
tal of 3145728 pixels separated on average by Spix = 6.87'. 
The maps have been cleaned of galactic foreground emis- 
sion using external templates, as explained in Gold et al. 
(2008). It is indeed necessary to mask out regions of strong 
foreground emission. Of the different masks provided from 
the LAMBDA Legacy Archive, which allow for the selective 
exclusion of portions of the sky at different flux levels, we 
chose the most conservative, KQ75 - essentially an exten- 
sion and improvement of the standard KPO cut in the 3-yr 
data release. This choice removes about 28% of the pixels. 
Point sources are masked based on a combination of exter- 
nal catalog data and WMAP-detected sources (Wright et al. 
2008). Note that point sources are the largest astrophysical 
contaminant to the temperature power spectrum. 

Our theory predictions (Section [2]) depend on the TT 
power spectrum, on the shape of the WMAP beam, and on 
the noise properties. We used the beam and window transfer 
functions from the five-year WMAP data, which have been 
significantly improved with the 5yr data release; models of 
the instrument gain and beam response are now accurate to 
better than ~ 0.6%, and errors in the five-year beam trans- 
fer functions are reduced by a factor of ~ 2 as compared 
to the three-year analysis (Hill et al. 2008). The WMAP5 
TT power spectrum we use comes from a weighted combi- 
nation of 153 individual cross-power spectra, and has been 
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improved by using a Gibbs-based maximum likelihood es- 
timate for I ^ 32 (Dunkley et al. 2008) and a pseudo Ci 
estimate for higher I (Nolta et al. 2008). The pixel-pixel co- 
variance matrix is essentially diagonal (Hinshaw et al. 2003, 
2008), so the rms noise in pixel p is a{p) — (tq/ A'obs(p), 
where the noise per observation, ctq, is provided for each 
coadded channel. This value for uo is obtained by averaging 
the values of the five possible year-by-year difference com- 
binations (Jarosik et al. 2007; Hinshaw et al. 2008). 

The data reduction process should be free of contami- 
nation induced by foregrounds, and should be insensitive to 
monopole and dipole moments. At A'aide = 512, the generic 
data pipeline involves the following main steps: 

(i) dipole and monopole removal outside KQ75; 

(ii) selection of pixels above a temperature threshold; 

(iii) HEALPix coordinate and pixel- noise assignment. 



3.2 Inhomogeneous Noise Properties 

The distribution of the rms noise values in a pixel at A'^side = 
512 is shown in MoUweide projection in Figure [T] and as 
histogram in Figure [2] for the WMAP5 Wl channel. The 
distribution is reasonably broad, indicating that the noise 
is inhomogeneous. This is a direct consequence of the fact 
that not all pixels were observed the same number of times. 
The evident symmetrical pattern present in Figure [T] is due 
to the scanning strategy of the WMAP satellite. 

In particular, equations p6|l and (I17|) show that the 
clustering of pixels separated by 9 depends on the joint dis- 
tribution oi p{(j\,(J2) given 9. Figure |3] shows an example of 
how the joint distribution p{(T\,02\9) varies as a function of 
separation 9 for the resolution Aside = 512. The four panels 
were constructed by randomly selecting 1000 pixels, and for 
each selected value of (ti, looking at the corresponding G2 
distributions at four different angular distances. Note that 
the scatter around the one-to-one line increases with increas- 
ing 9. This is not surprising; nearby pixels have similar Aobs 
whereas more widely separated ones do not (Figure [!}. Fig- 
ure|4]shows some examples of the corresponding conditional 
distributions p(o-2|cri, 0). When inserted into equations p6p 
and (|17p . such fits allow a prediction of how the clustering 
of pixels above a certain threshold temperature should de- 
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Figure 2. Distribution of the rms noise-per-pixel at Afgjdc = 512 
for the WMAP5 Wl coadded channel, after application of the 
KQ75 mask. 
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Figure 3. Joint distribution p{o"i, cF2\d) at Afgjdc = 512 resolution 
for four different choices of d as indicated in the panels, from the 
WMAP5 Wl coadded channel. 



pend on threshold if the signal is Gaussian and the noise is 
inhomogeneous. 



3.3 Distribution of Temperature 

The histogram in Figure [5] shows the distribution of temper- 
ature in the data, for the Wl coadded channel. For homo- 
geneous noise, the predicted distribution would be given by 
convolving the expected signal (in this case Gaussian with 
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e= 100', de= 3' 



0,12 




Figure 4. Examples of conditional probability distributions 
p((T2|o-i, 6»), when 6 = 100', for the WMAP5 Wl coadded channel. 
The panel shows results for N^^^^ = 512. Different sets of curves 
(spline fits) show the distribution of (T2 for tri =0.122, 0.142, 0.162 
mK, respectively. 



WMAP5 Wl coad, Kq75 mask 




-1 -0.5 0.5 1 

(5T [mK] 

Figure 5. Distribution of temperature when A^side = 512 for the 
WMAP5 Wl channel (histogram). Solid line shows the expected 
distribution (equation [TJ given the distribution of the noise. Dot- 
ted line is a Gaussian with the same rms as the data. 



variance determined by smoothing the best-fit power spec- 
trum on the scale of a pixel) with that of the noise (Gaus- 
sian with average rms ctn). For inhomogeneous noise, the 
predicted distribution (i.e. equation [1} is given by convolv- 
ing the expected Gaussian signal with the noise of rms an, 
and then weighting by the distribution of cr„ shown in Fig- 
ure [2l The dotted line in Figure [5] shows a Gaussian with 
the same rms as the data. This fit has shorter tails than 
the data. While this may be indicating that the signal is 
non-Gaussian (in fact, the sample skewness of the temper- 
ature distribution is —0.01216 mK^, and the sample excess 
Icurtosis 0.12698 mK*), the solid curve in Figure [S] shows the 
predicted distribution when inhomogeneous noise effects are 
included; it provides substantially better agreement with the 
measurements. 

Nevertheless, the one-point distribution itself cannot 
rule out a priori the presence of a non-Gaussian signal, for 
reasons made explicit in Fang & Pando (1997) ; hence it is not 
a suitable statistics to investigate departures from Gaussian- 
ity. This is also relevant for the detection of non-Gaussianity 
claimed by Jeong & Smoot (2007), where the one-point dis- 
tribution was used. We also note that their equation (9) is 
missing one integration over the distribution of the noise 
(see again our equation [1] for comparison). 

3.4 Pixel Number Density 

The number density of pixels above (below) a given tem- 
perature threshold is simply obtained by multiplying the 
fraction of pixels above (below) the threshold Dt (eq. |4]) 
by the total number density after masking (i.e. the total 
number of pixels after masking over the "available" or "un- 
masked" area). If the noise were homogeneous, then this 
prediction is rather straightforward: equation (Q simply re- 
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Wl 5yr coad, N = 512 
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Figure 6. Number density of hot and cold pixels at A'si^e = 512, 
for the WMAP5 Wl channel. Points are measurements from the 
data. Dotted line is the theoretical prediction for homogeneous 
noise, and solid line includes inhomogeneity. Different tempera- 
ture thresholds are considered, as explained in the main text. 

duces to erfc(i//\/2)/2, where u = Dt/o-u- In practice, due 
to the inhomogeneity of the noise, one needs to integrate Q 
numerically, using the distributions p{(t) measured from the 
data (Figure [2]). The solid lines in Figure |6] are the results of 
this integration, when accounting for inhomogeneous noise. 
Dotted hues are the theoretical predictions which assume 
homogeneous noise. Points in Figure |6] show measurements 
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Figure 7. Unweighted correlation functions of pixels above tlireshold calculated from the WMAP5 Wl temperature field, at A^sido = 512. 
A variety of pixel-thresholds are considered, as indicated in the panels. Dotted curves in the left panels show the predictions associated with 
Gaussian signal plus homogeneous noise, solid curves in the right panels show the case of inhomogeneous noise. Points are measurements 
of the clustering of hot and cold pixels, at corresponding temperature thresholds. Shaded areas are the Ic optimal variance errors. 



from the WMAP5 Wl channel, for pixels above or below a 
threshold, where \v\ = 0.5, 1.0, 1.5, 2.0, 2.5, respectively. Er- 
rorbars are from a Poisson analysis. While the solid lines 
slightly improve the fits, significant departures still remain, 
especially at higher thresholds. 

3.5 Two-point Statistics above Threshold 

We have measured the correlation function of pixels above 
and below a threshold temperature from the WMAP five- 
year coadded maps, for a range of thresholds. In each case, 
we measured the signal in two ways: one in which all pix- 
els are treated equally, and another in which the pixels are 
weighted by their temperature. We adopt here the standard 
estimator l~\-£^{9) — DD{9) / RR{8) for the correlation func- 
tion (Jensen & Szalay 1986), where DD{6) and RR{6) axe 
the number of data and random pairs, respectively. In partic- 
ular, the number of random pairs is computed by distribut- 
ing random points on a unit sphere, and then by applying the 
same procedures (i.e. KQ75 masking and HEALPix coordi- 
nate assignment) as for the data pipeline. For the weighted 
correlation estimates, the number of data pairs, DD{9), is 
simply replaced by its temperature weighted counterpart, 
WW{e). We also tested the Landy & Szalay (1993) and the 



Hamilton (1993) estimators, and found no significant differ- 
ences in the calculations. 

Figures[7]and[8]show an example of the results for pixels 
with \v\ = 0.25, 0.50, 0.75, 1.00, 1.25 and \v\ = 1.50, respec- 
tively. In all panels, symbols show the measurements of hot 
and cold pixels at the fixed threshold; errors are from a Pois- 
son analysis. Dotted curves in all the left panels show the 
results based on homogeneous noise, solid curves in all the 
right panels account for the fact that the noise is inhomo- 
geneous (equations 1121 and I13|l . Shaded areas are derived as 
explained in Appendix A (see also Hernandez-Monteagudo 
et al. 2004), and represent the la optimal variance errors. 

Accounting for inhomogeneity significantly improves 
the fit to the data. However, discrepancies between the ob- 
served correlation functions and the theoretical predictions 
still remain, especially for pixels at higher thresholds, or at 
small angular scales. In particular, we note the interesting 
fact that hot and cold pixels cluster differently within the 
same temperature thresholds. This feature is also present in 
all the other WMAP5 channels (see Appendix B), but we 
find that our theory is always in good agreement with the 
clustering of cold pixels in the Q and V frequency bands, and 
only in slight disagreement with the hot patch clustering, es- 
pecially for pixels at higher thresholds. Moreover, the fact 
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Figure 8. Same as Figure [7] but for the weighted correlation functions. 



that the WMAP5 data appear to be systematically stronger 
(weaker) at small (large) angular scales than the theoretical 
predictions, even at lower temperature thresholds (Figures[7] 
and|S]), is particularly true in the W channel (see Figure [B4l 
for comparison). We suggest some possible interpretations 
in the discussion section and in Appendix B, and present a 
more detailed investigation of this effect using realistic non- 
Gaussian mock simulations in a forthcoming paper. 

While the detected trend may be a unique signa- 
ture of primordial non-Gaussianity (i.e. a primordial non- 
Gaussianity would in fact enhance the clustering of the cold 
pixels and reduce that of the hot ones) , we study the effect of 
smoothing the map next, and provide some other plausible 
explanations for these anomalies in the final section. 

3.6 Effect of Smoothing 

We have investigated the effect of smoothing the WMAP 
five-year maps with a Gaussian beam. In particular, we have 
tried different smoothing scales, corresponding to a FWHM 
of 30', 45', 60' and 75', respectively. We have repeated the 
same analysis as for the unsmoothed maps, and character- 
ized the number density and the clustering statistics of pixels 
above threshold. 

Figure |9] summarizes the results for the number den- 
sity (Wl channel), after smoothing with different Gaussian 



beams. The number of pixels above or below the considered 
temperature threshold is reduced (a function of the smooth- 
ing scale), and measurements from the data and theoretical 
predictions appear now in agreement within about la at the 
threshold levels considered. 

Finally, we have characterized the clustering statistics 
and showed in Figure[Tn]the case of inhomogeneous noise, for 
an arbitrary choice of the threshold — 1.0). The Gaus- 
sian smoothing is included in our theory via equations ([24} 
and ((25}, in addition to the change of the noise variance dis- 
tribution. After smoothing the map with a Gaussian beam, 
the pixel noise may be no longer independent among dif- 
ferent pixels (i.e. {nin2) 7^ 0). However, accounting for this 
contribution (equation I10[) resulted in no significant differ- 
ence, even at small angular scales. In fact, the average vari- 
ance of the noise in the Wl channel is 0.0213 mK^ which 
is about 0.32% of the typical variance of the corresponding 
cosmological signal. Since the zero-level amplitude of the 
noise is already very low with respect to that of the signal, 
the mean amplitude of the correlated noise will be also very 
low. In addition to that, note that the WMAP mission has 
been designed to have uncorrelated pixel noise, and if such 
correlations arise when smoothing the maps, their amplitude 
is expected to drop quickly with the smoothing scale. 

When the smoothing scale is 75', we find that the ten- 
sion between data and theory is alleviated. 
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Figure 9. Number density of liot and cold pixels at N^\^^ = 512, for the smootlicd WMAP5 (Wl channel) maps. In each panel, points 
are measurements from data, dotted lines are the theoretical predictions for homogeneous noise, solid lines include inhomogeneity. From 
top to bottom and clockwise, the FWHM of the Gaussian beam is 75', 60', 45' and 30', respectively. 



4 DISCUSSION 

The observed distribution of temperatures in WMAP5 
pixels is slightly non-Gaussian (Figure [5]). This departure 
from Gaussianity is not unexpected if the intrinsic signal 
is Gaussian but the noise distribution is inhomogeneous 
white noise (Figures [T] and [Sjl. In fact, when accounting 
for inhomogeneous noise (Section [5}, we find that the 
one-point statistics is substantially in better agreement 
with the measurements. However, the predicted dependence 
of clustering on pixel height is not in good agreement with 
the data, even after accounting for the inhomogeneity of the 
noise (Figures [7l|8}. Although we found that an appropriate 
smoothing scale would be able to alleviate the discrepancy 
between data and measurements (Figures [9l and [TO)l . other 
possible reasons for these discrepancies are: 

Masking procedure. In our analysis we adopted the 
KQ75 mask, which allows for selective exclusions of bright 
portions of the sky (28% pixel cut). Although this is rather 
a drastic cut and the mask was significantly improved with 



the WMAP five-year data release (Hill et al. 2008) , it may 
not completely account for all the Galactic foreground ef- 
fects. However, cutting out all the Galactic plane (|6| ^ 30° 
strip) and repeating our analysis resulted in no significant 
change. Edge effects due to pixels which lie very close to the 
mask could also affect our analysis. 

Cold spot contamination. Inoue & Silk (2007) suggest that 
the presence of low density regions in the southern Galac- 
tic cap could account for our anomalies, but restricting our 
analysis to b ^ 30°, the northern Galactic cap, resulted in 
no significant change. 

Contamination by point sources. Lopez-Caniego et al. (2007) 
detected 98 new sources (i.e. 26%) not present in the WMAP 
three-year catalog. At higher frequencies WMAP estimates 
neglect the deviations of the point spread function from a 
Gaussian shape, and carry out a blind search for sources. 
With the WMAP five-year data release, the mask for point 
sources has been significantly improved (Wright et al. 2008), 
and other studies (Gonzalez-Nuevo et al. 2008; Chen & 
Wright 2008) updated the point source catalog. However, 
Massardi et al. (2008) recently provided a new version of 
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the catalog based on WMAP5 data and found new point 
sources (484 sources detected), carrying a complementary 
blind and non-blind approach. This may be important to 
our analysis, since even a very low level contamination in 
the maps can produce spurious non-Gaussianities. In fact, 
we find some discrepancies at small angular scales (around 
20') especially in the W band, and their origin may be asso- 
ciated with extragalactic dust emission (one type of point- 
source contamination), which peaks at high frequencies. In 
a forthcoming study, we will be addressing the effect of con- 
tamination induced by point sources in more depth, using 
an updated version of the source catalog and mock simula- 
tions. 

Foreground subtraction contamination. Uncertainties in the 
external Foreground Template Model used for the fore- 
ground subtractions (Gold et al. 2008) may introduce 
anomalies at the percentage level. The template itself has 
noise, which may be correlated at small separations. This is 
a delicate issue, since galactic foregrounds are non-Gaussian 
and anisotropic, and even low level contamination in the 
maps can produce detectable non-Gaussianities (Park et al. 
2002; Naselsky et al. 2005; Kim, Naselsky & Christensen 
2008). On the other side, it is also worth noticing that re- 
cently Vio & Andreani (2008) showed that the benefits of 
using more sophisticated methods for foreground cleaning, 
such as the Harmonic Internal Linear Combination, are over- 
estimated. We are also addressing the foreground subtrac- 
tion systematics in a forthcoming study. 
WMAP beam, window function, absolute calibration and 
pixel-noise uncertainties. Window functions were computed 
from the symmetrized beam profiles following the Hermite 
method in Page et al. (2003). A typical WMAP window 



function has an uncertainty of 2 — 3% (Hill et al. 2008), and 
these uncertainties add in quadrature in the cosmological 
analysis. Absolute calibration uncertainties in the five-year 
WMAP data are estimated to be 0.5% (Hinshaw et al. 2008). 
Recent studies (Groeneboom et al. 2009; Kathrine Wehus et 
al. 2009) have also pointed out that the noise levels of these 
maps are underestimated, and that there are some problems 
with the standard WMAP transfer functions as well. It may 
be that these facts play an important rule in our study (see 
also Colombo et al. 2008). 

Secondary anisotropics and phase transitions in the early 
universe. Spurious non-Gaussianities could arise from sec- 
ondary anisotropics, such as gravitational lensing, Sunyaev- 
Zel'dovich effect or Sachs- Wolfe effects (Babich & Pierpaoli 
2008; Carbone et al. 2008). Phase transitions in the early 
universe may also introduce a new source of non-Gaussianity 
(Silvestri & Trodden 2008). All these effects may be difficult 
to disentangle from a pure primordial non-Gaussian signal, 
even after a clear detection of primordial non-Gaussianity. 
Real non-Gaussian signatures at small scales. A better un- 
derstanding of all the previous points is necessary before we 
can claim that the signal we see is due to primordial non- 
Gaussianity. However, we found an interesting difference in 
the clustering of hot and cold pixels within the same temper- 
ature threshold level; this fact may be a unique signature of 
primordial non-Gaussianity. We present a detailed investi- 
gation of this feature in a forthcoming study. We finally note 
that, as this work was being refereed, Hon et al. (2009) pre- 
sented a frequentist analysis of the correlation functions of 
the local extrema, and also found inconsistency with Gaus- 
sian simulations plus differences in the clustering of hot and 
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cold peaks - although comparison with their work is not 
direct since they considered bigger angular scales. 

While seeking for primordial non-Gaussianity is at the 
moment a new frontier in cosmology, ongoing efforts are 
currently devoted to the characterization of non-Gaussian 
confusion effects, to reliable theoretical predictions of non- 
Gaussianity from models (see for example Boyle & Stein- 
hardt 2008; Fergusson & Shellard 2008; Munshi & Heav- 
ens 2009), till the extraction of information from data (i.e. 
Raeth et al. 2008; Gong et al. 2009) or the search for obser- 
vational signatures of primordial non-Gaussianity imprinted 
in the large-scale structure of the universe. Even a small de- 
gree of primordial non-Gaussianity can be a crucial probe 
of the inflationary dynamics or alternative universe models, 
hence studies of non-Gaussianity may eventually become a 
powerful and solid probe of ultra-high energy physics and 
inflation. 

We note that our model for the effects of inhomogeneous 
noise may be useful in other studies (see for instance Yu & 
Lu 2008). Extending our formalism for inhomogeneous noise 
to peak rather than pixel statistics (for example Heavens & 
Sheth 1999; Heavens & Gupta 2001) is more complicated; 
this will complement numerical Monte-Carlo analyses of this 
problem (Larson & Wandelt 2004; Tojeiro et al. 2006; Ayaita 
et al. 2009; Hou et al. 2009), and is the subject of work in 
progress. 
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APPENDIX A: ACCOUNTING FOR FINITE 
BINSIZE EFFECTS 

We estimate the uncertainties in the correlation functions 
above threshold from the optimal variance limit, contain- 
ing cosmic variance, instrumental noise, and finite binsize 
effects. 

In more detail, the ultimate accuracy with which the 
CMB power spectrum can be determined at each / is given 
by (Knox 1995): 



'^ov(0 



(2; + l)/,ky 



N 14/WMAP 



ACi (Al) 



where is the number of pixels, ctn the average pixel noise, 
^WMAP ^j^g WMAP window function, and /sky the fraction 
of the sky covered by the experiment. The uncertainty in 
the angular correlation function for narrow bins in 9 is then 
(Hernandez-Monteagudo et al. 2004): 
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where Ci = l{l + l)Ci/27r and ig the additional 

smoothing due to finite pixel size, mask infiuence and an 
optional Gaussian beam. 

However, in practice C{9) is not measured at a point, 
but smeared out over a region of size A9. In the limit of large 
I and smaU 9, Pi{cos9) ~ Jo[{l + 1/2)9] (Bond & Efstathiou 
1987). If the bin size is not infinitesimal, one must replace: 



Dotted lines in the figure are Gaussians with the same rms 
as data, solid curves show the predicted distributions when 
inhomogeneous noise effects are included. Departures from 
the Gaussian fits are more significant in the W frequency 
range. Figure lB3l shows the number density of pixels above 
(below) different temperature thresholds. Dotted lines are 
theoretical predictions with homogeneous noise, solid lines 
include inhomogeneity. Points are measurements from all the 
WMAP5 DA's. Finally, Figure lB4l highlights some examples 
of the unweighted correlation functions from the WMAP5 
pixel-pixel temperature fields, at Agidc ~ 512. Two pixel- 
thresholds are considered = 0.25 and \u\ — 0.75), as 
indicated in the panels. Solid curves show the predictions 
associated with Gaussian signal plus inhomogeneous noise. 
Points are measurements of the clustering of hot and cold 
pixels at corresponding temperature thresholds, and shaded 
areas are the la optimal variance errors. Essentially, we find 
that the difference in clustering between hot and cold pixels 
is still present in all the other DA's. However, our theory is 
always in good agreement with the clustering of cold pix- 
els, and in slight disagreement with that of the hot ones, 
both in the Q and in the V channels. Therefore the 94 GHz 
frequency (W band) seems to be the most discrepant one, 
with respect to our theoretical predictions. In particular, we 
checked our small-scale theory expectations against WMAP 
mock measurements from simulations with identical beam 
and noise properties, and found full consistency. Hence, we 
suggest that those small-scale discrepancies may be due to 
extragalactic dust emission (one type of point-source con- 
tamination), which has a peak at high frequencies and is 
not accounted in our analysis. 
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where t-mi-n ~ 9~A9/2 and tmax ~ 9 + A9/2. This expression 
should then substitute Pi (cos 61) ~ Jo[{l + 1/2)9] in (|X2|| . to 
give o"c^g(9) = AC AO (9). The uncertainties in the correla- 
tion function above threshold are finally derived from: 
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APPENDIX B: RESULTS FOR ALL THE 
OTHER WMAP FIVE-YEAR DIFFERENCING 
ASSEMBLIES 

We have inspected all the WMAP five-year channels (W 
band at 94 GHz; V band at 61 GHz; Q band at 41 GHz), 
and performed the same analysis described in the main text 
for each single differencing assembly (DA). For this study, it 
is ideal to consider individual DA's rather than their com- 
binations, as WMAP beams and noise properties are well- 
defined within each single DA. Figure IBll shows the dis- 
tributions of the rms noise values for all the DA's consid- 
ered; note the characteristic bi-modality of these histograms, 
present in all the channels. Figure IB2I shows the distribu- 
tions of temperatures in the data at different frequencies. 
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Figure Bl. Distributions of the rms noise-per-pixel at A^sidc = 512 for the five-year coadded WMAP data, after apphcation of the KQ75 
mask. Prom top to bottom and from left to right, the various differencing assemblies are respectively Ql and Q2 (the 41 GHz channel), 
VI and V2 (the 61 GHz channel), and W2, W3, W4 (the 94 GHz channel). 
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Figure B2. Distribution of temperatures when A^side = 512 (histograms). Solid lines in all panels show the expected distributions 
(equation [TJ given the corresponding distributions of the noise. Dotted lines are Gaussians, with the same rms as the data. Prom top to 
bottom and from left to right, the various DA's are respectively Ql, Q2, VI, V2, W2, W3, W4. 
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Figure B3. Number densities of hot and cold pixels at N^^^^ = 512, for the WMAP5 individual channels. Points are measurements from 
the data. Dotted lines in all the panels are theoretical predictions for homogeneous noise, solid lines include inhomogeneity. From top to 
bottom and from left to right, the various DA's are respectively Ql, Q2, VI, V2, W2, W3, W4. 
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Figure B4. Examples of unweighted correlation functions from the WMAP5 pixel-pixel temperature field, at A^sidc = 512. Two pixel- 
thresholds are considered = 0.25 and = 0.75), as indicated in the panels. Solid curves in all panels show the predictions associated 
with Gaussian signal plus inhomogeneous noise. Points are measurements of the clustering of hot and cold pixels, at corresponding 
temperature thresholds. Shaded areas are the Icr optimal variance errors. From top to bottom and from left to right, the various DA's 
are respectively Ql, Q2, VI, V2, W2, W3, W4. 



